Generic behaviour of nonlinear sound waves near the surface of a star: smooth 

solutions 
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We are interested in the generic behaviour of nonlinear sound waves as they approach the surface 
of a star, here assumed to have the polytropic equation of state P = Kp v . Restricting to spherical 
symmetry, and considering only the region near the surface, we generalise the methods of Carrier 
and Greenspan (1958) for the shallow water equations on a sloping beach to this problem. We give 
a semi-quantitative criterion for a shock to form near the surface during the evolution of generic 
initial data with support away from the surface. We show that in smooth solutions the velocity and 
the square of the sound speed remain regular functions of Eulerian radius at the surface. 



I. INTRODUCTION 

In numerical simulations of neutron stars in general 
relativity, the matter is often modelled as a perfect fluid. 
The simplest equation of state usually considered is the 
ideal gas equation of state 



P = (r - l)ep, 



(1) 



where P is the pressure, p the rest mass density and e 
the internal energy per rest mass. The polytropic index 
r=l + l/n>lisa constant. 

If the entropy per rest mass is everywhere the same, 
the ideal gas equation of state reduces to the polytropic 
equation of state 



P = Kp L 



(2) 



where K is another constant depending on the entropy 
per rest mass. If the initial data are isentropic, the solu- 
tion remains isentropic until a shock forms. For the poly- 
tropic equation of state with n > 0, spherically symmet- 
ric self-gravitating solutions with a regular centre (stars) 
have a surface, characterised by P = p = at finite ra- 
dius r = r*, where p ~ (r* — r) n near the surface. 

Standard numerical methods for evolving stars fail at 
the surface because division by zero density occurs and 
the speed of sound goes to zero. For smooth solutions 
in spherical symmetry, this can be avoided by using La- 
grangian coordinates, but in 3-dimcnsional (3D) simula- 
tions with high-resolution shock capturing (HRSC) meth- 
ods, the standard practice is to match the star to a thin 
"atmosphere" , which is then artificially kept from ac- 
creting onto it. This method is likely to give qualita- 
tively wrong results, as the wave structure of the Rie- 
mann problem that underlies HRSC methods is different 
if the right state is vacuum. 

The failure of the numerical methods is related to the 
physical fact that the perfect fluid approximation must 
break down at the surface. This approximation includes 
the assumption that small fluid elements are in thermal 
equilibrium on dynamical timescales, but as the density 
goes to zero, the thermal timescale diverges while the 
fluid dynamical timescales are still determined by waves 



in the interior and remain finite. In reality, some kind of 
plasma physics approximation applies. 

The premise of this paper is that a mathematically 
correct numerical implementation of the perfect fluid as- 
sumption is more correct than the use of an unphysical 
atmosphere, which at best introduces physically unmo- 
tivated approximations and at worst does not even have 
a continuum limit. In this paper we provide two math- 
ematical results that should be useful in achieving this 
goal. We begin here with smooth solutions and leave 
shocks for later work. 

Our preliminary question is whether smooth initial 
data representing an outgoing wave with compact sup- 
port form a shock as the wave approaches the surface. 
That a shock forms is suggested by the fact that the 
sound speed goes to zero at the surface with c s ~ \/r* — r 
(independently of the polytropic index n), so that any 
outgoing wave steepens. Sperhake [1] has investigated 
this numerically in general relativity in spherical sym- 
metry and concludes that small amplitude waves do not 
shock but large amplitude waves do. In the Newtonian 
case in spherical symmetry this had already been proved 
by Pelinovsky and Petrukhin [3] . We improve on this re- 
sult by deriving a semiquantitative criterion for a sound 
wave to remain regular as it approaches the surface. 

Our main question is what kinematic boundary condi- 
tions can be used in a numerical simulation to represent 
the free boundary at the surface of the star. This has 
been addressed in general relativity by Sperhake [1] for 
nonlinear spherical perturbations, using Lagrangian co- 
ordinates, and by Passamonti [2] for linear non-spherical 
perturbations. Here we consider the nonlinear case in 
Eulerian coordinates. 

To answer both questions we use the mathematical 
methods of a classic paper by Carrier and Greenspan 
[4] concerning the shallow water equations on a sloping 
beach. We begin by reviewing their results and extending 
them from the shallow water case n — 1 to the general 
polytropic case n > 0. 
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II. MATHEMATICAL SETUP 

For simplicity we assume spherical symmetry. Near 
the surface of the star, gravity is typically weak. Fur- 
thermore, the formation of shocks does not require large 
fluid velocities. This suggests that Newtonian physics 
should be a good approximation for what we want to 
investigate. On a sufficiently small scale the spherical 
symmetry of the star reduces to planar symmetry, and 
the Newtonian gravitational acceleration g is dominated 
by the interior of the star, and can be approximated as 
constant in space and time. Finally, for smooth solutions, 
sufficiently close to the surface, the entropy gradient can 
be neglected compared to the density gradient in deter- 
mining the pressure gradient. We can therefore approx- 
imate the ideal gas as isentropic, with equation of state 
(2). (This last approximation would not hold if a shock 
reached the surface.) 

In the "radial" spatial coordinate x and time t, with v 
the Eulerian fluid velocity in the x direction, the Euler 
and conservation equations are 

v t + vv x + TKp 1 ^- 1 p x = -g, (3) 
Pt + vPx + pv x = 0. (4) 

Here p = defines a free boundary x — x*(t). Within 
the approximation of planar symmetry, x has an infinite 
range, with x < x*(t) representing the interior of the 
star. 

It is useful to replace the dependent variable p with the 
sound speed c given by c 2 = dP/dp = YKp x l n to obtain 

v t + vv x + 2ncc x = -g, (5) 
I 

Ct + vc x + -^cv x = 0. (6) 

For n = 1, these equations are identical with the shal- 
low water equations restricted to planar symmetry on 
a uniformly sloping beach, with x and v the horizontal 
position and velocity, p ~ c 2 the height of the water, 
g the effective horizontal gravitational acceleration, and 
x = x* (t) the instantaneous shoreline [5] . 

The unique static solution of (5,6) is 



v = 0, c = 




and hence p <~ (— x) n , where we have fixed a translation 
invariance by locating the surface at x = 0. 



(considering the shallow water case n = 1) suggested a 
hodograph transform from independent variables t and x 
to independent variables A and a given by 

A = v + gt, (9) 
a = 2nc. (10) 

(These definitions differ from [4] by a factor of 2.) 
The resulting transformation of partial derivatives is 




where 

A = t x x a - x x t a . (12) 
In particular, we have 

(£::)=*-'(-*; T> < i3 > 

Clearly, the transformation is regular if and only if A ^ 
0,±oo. 

Substituting (9-10) and (13) into (8), we obtain 

x a - (A - gt) t a + (^) t x = 0, (14) 
xx+(^)t a -{X-gt)t x = 0. (15) 

This PDE system is not yet linear because of the appear- 
ance of gt in the coefficients of t a and t\. However, from 
the two nonlinear first-order PDEs (14-15) one can derive 
a linear second-order PDE for i(A, er), namely 

2n + 1 l n\ 

t\\ = t aa H t(j . (16) 

a 

Trivially, A taken as a function of A and a obeys the 
same PDE as t, and by adding the two we obtain the 
autonomous linear wave equation 

2n + l 

v\x = iw H v a (17) 

(7 

for u(A, a). This is the key equation of this paper. 

The problem has now been cast into linear form, and 
the free boundary x = x* (t) has been mapped to the co- 
ordinate line cr = 0, with a > representing the interior 
of the star. 



III. HODOGRAPH TRANSFORM IV - A CRITERION FOR SHOCK FORMATION 

The problem can be written as From ( 13 ) with ( 9 > 10 ) wc find 

[d t + (v±c)d x ]{v + gt±2nc)=0, (8) Vx = _L V(Tj (18) 

and so admits the Ricmann invariants (v + gt)±2nc with 1 

characteristic speeds v ± c. Carrier and Greenspan [4] c ' x ~ 2ngA ' 
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and so a shock forms from regular initial data as and only 
if A — > 0. Using (9-10,14-15), the Jacobian A defined by 
(12) can be expressed in terms of v alone as 



A 



2ng 2 



(l-v x f-vl 



(20) 



We see that the wave does not form a shock if the first 
derivatives of v in a solution of (17) remain sufficiently 
small, so that A remains negative. Such solutions are 
easily obtained by rescaling the amplitude of any given 
solution. 

We shall now consider small smooth initial data for 
(5,6) on the curve t = 0, x < 0. These correspond to 
Cauchy data for (17) on the curve given by A = Ao(ct), 
a > 0. We require these data to obey 



(1-«a)' 



> 



(21) 



for all p > on A = A (cr). This criterion is necessary 
for the existence of the equivalence between (5,6) and 
(17), and implies that there is no shock present in the 
initial data. We then formally evolve the data to A > 
A (er) using (17). Setting aside the boundary at a = 0, 
which we consider later, this solution exists because (17) 
is linear. However, if at any point in A > Ao(cr) the 
condition (21) is violated, the wave has developed a shock 
at some t > 0, and the solution of (17) does not have 
physical meaning for larger values of t. 

In order to translate initial data in coordinates (x, t) 
to (a, A) , we consider smooth data with compact support 
away from the boundary and which are sufficiently weak 
(in the sense of close to the static star solution) that ini- 
tially the solution can be approximated by a solution of 
the linearisation of (5,6) around the static star solution. 
We then evolve these data using (17), and so do not re- 
quire them to remain small. We use (21) in this solution 
as the necessary and sufficient criterion for the absence 
of shocks. 

Linearising (5,6) about the static solution (7), we ob- 
tain 



SVf. 



("?) 



Sv x 



1 



-Sv a 



(22) 



(We have written Sv instead of v to stress that this is only 
an approximation valid for small v.) The same equation 
can be obtained from (17) by the substitutions 



A = gt, 



2^—gnx. 



(23) 



This gives us a simple approximate relation between ini- 
tial data for the linearisation of (5,6), and initial data for 
(17) (which is linear but contains the nonlinear dynam- 
ics). 

A formal d'Alembert solution of (17) is h 



"^(Aia) 



(24) 



where is free data and 



fk+l 



(k 



2{k+l) 



I 



ft- 



(25) 



A few remarks will put this result into context: This 
series converges at most in the sense of an asymptotic se- 
ries as o — > oo, and clearly diverges for sufficiently small 
a. Another formal d'Alembert solution exists which has 
ascending powers of a, but it does not interest us here. In 
the special case n = 1/2, either series reduces to the well- 
known d'Alembert solution of the spherical wave equa- 
tion in 3 dimensions, while for n = —1/2 we obtain the 
d'Alembert solution of the 1-dimcnsional wave equation. 

Consider now an isolated wave packet approaching the 
surface with initial position cro, width o\ <C o~o and am- 
plitude vq, so that \v\\ ~ \v a \ ~ vq/<ji initially. In this 
regime, we can approximate 



v(\,o-)~a- n -*f+(\ + (T). 



(26) 



The derivatives of v take their largest values when the 
wave packet turns around close to the surface. From 
the scaling properties of solutions of (17), this must hap- 
pen at <7 ~ oi, at which point its amplitude will be 
w o( cr i/ cr o)~ Il ~ 1/ ' 2 in the approximation (26). Evaluating 
(21) at that point, we obtain a criterion for the wave 
never to form a shock, which is 



vo < / cti 



o~ 1 



(TO 



(27) 



Finally, expressing oo and o\ in terms of the initial Eu- 
lerian position xo and length scale x\ of the wave packet 
by using (23), we obtain the regularity criterion 




n+(3/2) 



(28) 



In these estimates we neglect an unknown 0(1) factor 
depending on the precise shape of the wave packet. 

Although we have worked in the approximation of pla- 
nar symmetry and constant g, it is useful to express the 
parameter g in terms of = \/2gr\, which is the es- 
cape velocity at the surface of a spherical star, where r* 
is its radius and g is the gravitational acceleration at its 
surface. We can then rewrite the estimate (28) as 



'"o 



< 



xi 



-(3/2) 



Fol 



(29) 



m a neu- 

-|4„ 



± fe=0 



A numerical example will illustrate this: 
tron star modelled as a polytrope with r* ~ 10 4 m, 
v* ~ 10 8 to/s and n = 1, a sound wave of wavelength 
~ lm deep in the interior (xq ~ — r*) must have an 
amplitude of vq < 10~ 2 m/s to remain regular. 
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V. GENERIC BEHAVIOUR AT THE FREE 
BOUNDARY 

The surface of the star is a free boundary characterised 
by the kinematic boundary conditions 

P{x.(t),t) = 0, (30) 
^ = v(x.(t),t). (31) 

These conditions are straightforward to implement in La- 
grangian coordinates, but in 3D HRSC simulations we 
need their equivalent in Eulcrian coordinates. For so- 
lutions which remain smooth, we obtain these by going 
through the hodograph transformation. 

The general solution of Eq. (17) can be written as a 
linear superposition of solutions of the form 

v(\, a) = e wX o- n J± n (uj(j). (32) 

As J n (a) is a n times a power series in positive even pow- 
ers of a, the solution using J n is an even regular function 
of a, while the solution using J_„ diverges as a~ 2n as 
a — > 0. The regular solution can be selected by imposing 
the boundary condition 

v a = at cr = 0, (33) 

which together with (17) makes a well-posed linear 
initial-boundary value problem. Clearly this condition 
is the required kinematic boundary condition for smooth 
solutions. 

We now translate this back into the Eulerian variables 
c(x,t) and v(x,t). Assuming the wave does not form a 
shock, the square bracket in (20) is strictly positive, and 
so A <~ a at the boundary. Substituting A ~ a into (19) 
gives 

ac x ~ (1 - v x ) (34) 

at the boundary. The right-hand side is even in a because 
v is even in a by the assumption of regularity. It follows, 
using (10), that (c 2 ) x is a regular function of c 2 , and 
hence c 2 is a regular function of x. 
Substituting A ~ a into (18) gives 

v x ~ <y~ x v„ (35) 

at the boundary. The right-hand side is again even in 
a. Hence v x is an even function of c 2 and so, using our 
previous result, it is a regular function of a;. It follows 
that v is a regular function of x. 

It is clear that the A or t dependence does not affect 
these results in the limit x — > a;*(i) or a — > 0. We have 



therefore shown that as long as the solution remains reg- 
ular, c 2 and v are regular functions of x and t at the 
surface. This is the desired kinematic free boundary con- 
dition. In particular, c 2 ~ x*(i) — x at the moving surface 
of regular solutions, as in the static case. 

Note that p <~ (c 2 )", so p is a regular function of x 
only if n is an integer. Note also that in general v and c 2 
are neither even nor odd ini-i,. 

VI. DISCUSSION 

Building on the earlier work [3, 4], we have given vari- 
ous forms of an upper limit on the amplitude of nonlinear 
sound waves if they are to avoid forming a shock. This 
tells us in which physical regime a simple (non-shock cap- 
turing) numerical method will be valid because shocks do 
not occur. It may also be of direct astrophysical interest. 

For solutions which remain regular as they are reflected 
at the free boundary, we have shown that the usual free 
boundary condition is equivalent to v and c 2 being reg- 
ular functions of x and t. This suggests an alternative 
numerical treatment of the stellar surface which does not 
require an unphysical atmosphere. 

Our results were derived within the approximations of 
Newtonian physics, a constant gravitational field, a poly- 
tropic equation of state and planar symmetry (as the 
limit of spherical symmetry near the surface). As dis- 
cussed in the introduction, these arc all natural approx- 
imations to make, except for spherical symmetry. How- 
ever, applying geometric optics to the linearised sound 
wave equation for the pressure perturbation SP, 

SPtt = (-^) (sp xx + ^SP X + SP yy + SP zz ^j , 

. (36) 

we find that its sound rays, without loss of generality re- 
stricted to the xy plane, are given by y(x) = a+ln(l-|-6a; 2 ) 
for constants a and b, and so in the geometric optics ap- 
proximation sound waves approaching the surface x = 
at any angle are refracted towards lower sound speed un- 
til they reach the surface at right angles. This provides 
some justification for the assumption that our results will 
also be qualitatively correct beyond the restriction to 
spherical (planar) symmetry. 
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